Descriptive microbiome stats and cleaner code for dissertation / publication.

1 Load libraries and directories

library(magrittr)
library(phyloseq)
library(fs)
library(tidyverse)
library(vegan)
library(rstatix)
library(ggpubr)
library(microViz)
library(pals)
library(RColorBrewer)
library(colorBlindness)
library(microbiome)
library(ggrepel)
library(DT)
library(plotly)
library(ape)
library(picante)
library(phytools)
library(DESeq2)
library(parallel)
library(doBy)
library(ragg)
library(cowplot)
library(ggtext)
rm(list=ls()) # clear environment 

git.dir = file.path("/hpc/group/dmc/Eva/GLP_KO_Microbiome_final") # CHANGE ME 
ps.rds <- file.path(git.dir, "ps.GLP1.rds")
ps <- read_rds(ps.rds)
ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3304 taxa and 234 samples ]
## sample_data() Sample Data:       [ 234 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq()      DNAStringSet:      [ 3304 reference sequences ]
fig.dir = file.path(git.dir, "figures")

map.file = file.path(git.dir, "Metadata_GLP1.csv")
meta.df = read_csv(map.file, show_col_types = FALSE)

Note that this phyloseq object is already filtered such that only samples with reads > 1000 are available. See here:

sample_sums(ps) %>% as.data.frame() %>% 
  dplyr::rename(Reads = ".") %>% 
  arrange(Reads) 

The names are hard to interpret because they are the SRA accession numbers. Let’s change the names in the phyloseq object.

phyloseq::sample_names(ps) <- phyloseq::sample_data(ps)$Label

1.1 Change metadata to make variables factors

meta.df$Genotype <- factor(meta.df$Genotype, levels = c("WT", "KO"))
meta.df$Sex <- factor(meta.df$Sex, levels = c("Female", "Male"))
meta.df$Diet <- factor(meta.df$Diet, levels = c("Normal_Chow", "High_Fat_Diet"))
meta.df$Intranasal_Treatment <- factor(meta.df$Intranasal_Treatment, levels = c("PBS", "HDM"))
meta.df$Surgery <- factor(meta.df$Surgery, levels = c("None", "Sham", "VSG"))
meta.df$Group <- factor(meta.df$Group, levels = c("NA", "Control", "1", "2"))
meta.df %>% dplyr::filter(Label %in% rownames(phyloseq::sample_data(ps))) %>%
  mutate(Label_graph = Label) %>% 
  column_to_rownames("Label") -> meta.ps

ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3304 taxa and 234 samples ]
## sample_data() Sample Data:       [ 234 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq()      DNAStringSet:      [ 3304 reference sequences ]
sample_data(ps) <- meta.ps
ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3304 taxa and 234 samples ]
## sample_data() Sample Data:       [ 234 samples by 40 sample variables ]
## tax_table()   Taxonomy Table:    [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq()      DNAStringSet:      [ 3304 reference sequences ]

Genotype, Sex, Diet, Intranasal_Treatment, and Surgery are appropriately factors now.

1.2 Subset phyloseq objects

# Only include samples 
subset_samples(ps, Type == "True Sample") -> ps_sam

# remove samples with potential labeling issues 
ps_sam %>% subset_samples(row.names(ps_sam@sam_data) != "294" & row.names(ps_sam@sam_data) != "30") -> ps_sam

# Subset by sample type
ps_sam %>% subset_samples(Sample_type == "Lung Tissue") -> ps.LT.asv
ps_sam %>% subset_samples(Sample_type == "Fecal") -> ps.F.asv

ps_sam
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3304 taxa and 221 samples ]
## sample_data() Sample Data:       [ 221 samples by 40 sample variables ]
## tax_table()   Taxonomy Table:    [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq()      DNAStringSet:      [ 3304 reference sequences ]
ps.LT.asv
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3304 taxa and 44 samples ]
## sample_data() Sample Data:       [ 44 samples by 40 sample variables ]
## tax_table()   Taxonomy Table:    [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq()      DNAStringSet:      [ 3304 reference sequences ]
ps.F.asv
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3304 taxa and 177 samples ]
## sample_data() Sample Data:       [ 177 samples by 40 sample variables ]
## tax_table()   Taxonomy Table:    [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq()      DNAStringSet:      [ 3304 reference sequences ]
# Relative abundance
ps_sam.ts <- transform_sample_counts(ps_sam, function(x) x/sum(x))
ps.LT.asv.ts <- transform_sample_counts(ps.LT.asv, function(x) x/sum(x))
ps.F.asv.ts <- transform_sample_counts(ps.F.asv, function(x) x/sum(x))

1.3 Clean-up Metadata

Let’s subset the metadata to those that are in the ps_sam object.

meta.df %>%  
  filter(meta.df$Label %in% rownames(ps_sam@sam_data)) -> meta.sam

meta.sam %>% filter(is.na(Genotype) == TRUE & is.na(Intranasal_Treatment) == TRUE & is.na(Surgery) == TRUE)

2 Rarefaction curves

rarecurve(as.data.frame(otu_table(ps_sam)), step = 100, cex = 0.5, tidy = TRUE) -> rare_ps #put in ggplot2 friendly dataframe
## Warning in rarecurve(as.data.frame(otu_table(ps_sam)), step = 100, cex = 0.5, :
## most observed count data have counts 1, but smallest count is 2
meta.sam %>% dplyr::select(Label, Sample_type, Mouse, Week, DNA_ng_ul_raw) %>%
  dplyr::rename(Site = Label) %>% 
  right_join(rare_ps, by = "Site") -> rare_ps

ggplot(data=rare_ps, aes(x=Sample, y=Species, group=Site)) + 
  geom_line() + theme_bw() + labs(x = "Number of Reads", y= "Species Richness") + 
  facet_wrap(~Sample_type, scales = "free")  -> p
p

ggplotly(p, tooltip = "Site")

3 Relative abundance plots

ps_sam %>% tax_glom("Phylum") %>% 
  transform_sample_counts(function(x) x/sum(x)) -> ps.phylum.ts

taxa_names(ps.phylum.ts) <- tax_table(ps.phylum.ts)[, 'Phylum']
ps.phylum.ts
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 20 taxa and 221 samples ]
## sample_data() Sample Data:       [ 221 samples by 40 sample variables ]
## tax_table()   Taxonomy Table:    [ 20 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 20 tips and 19 internal nodes ]
## refseq()      DNAStringSet:      [ 20 reference sequences ]
plot_bar(ps.phylum.ts,x = "Label_graph", fill="Phylum") + 
  geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +  
  theme_bw()  +   
  labs(y="Relative Abundance", x = "Sample") +
  coord_cartesian(ylim = c(0,1), expand=0) + 
  scale_color_manual(values=SteppedSequential5Steps) + 
  scale_fill_manual(values=SteppedSequential5Steps) +
  theme(axis.text.x=element_blank()) +
  facet_wrap(~Sample_type, scales = "free") -> p

p

ps_melt(ps.phylum.ts)  -> phylum.ts.melt
phylum.ts.melt %>%
  filter(Sample_type == "Lung Tissue") %>% 
  doBy::summary_by(Abundance ~ OTU, FUN = median) %>%
  arrange(desc(Abundance.median))
phylum.ts.melt %>%
  filter(Sample_type == "Fecal") %>% 
  doBy::summary_by(Abundance ~ OTU, FUN = median) %>%
  arrange(desc(Abundance.median))
ps.phylum.ts %>% 
  subset_samples(Sample_type == "Fecal") %>% 
  ps_arrange(desc(Firmicutes), desc(Bacteroidota), .target = "otu_table") -> ps.fecal.phylum.ts

sample_data(ps.fecal.phylum.ts)$Sample_Order <- c(1:nrow(sample_data(ps.fecal.phylum.ts)))
ps.fecal.phylum.ts
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 20 taxa and 177 samples ]
## sample_data() Sample Data:       [ 177 samples by 41 sample variables ]
## tax_table()   Taxonomy Table:    [ 20 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 20 tips and 19 internal nodes ]
## refseq()      DNAStringSet:      [ 20 reference sequences ]
ps.phylum.ts %>% 
  subset_samples(Sample_type == "Lung Tissue") %>% 
  ps_arrange(desc(Proteobacteria), desc(Actinobacteriota), .target = "otu_table") -> ps.lt.phylum.ts

sample_data(ps.lt.phylum.ts)$Sample_Order <- c(1:nrow(sample_data(ps.lt.phylum.ts)))
ps.lt.phylum.ts
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 20 taxa and 44 samples ]
## sample_data() Sample Data:       [ 44 samples by 41 sample variables ]
## tax_table()   Taxonomy Table:    [ 20 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 20 tips and 19 internal nodes ]
## refseq()      DNAStringSet:      [ 20 reference sequences ]
ps_melt(ps.fecal.phylum.ts) %>%
  ggplot(aes(fill=Phylum, y=Abundance, x=Sample_Order)) + 
    geom_bar(position="fill", stat="identity", width = 1) +
  theme_bw()  +   
  labs(y="Relative Abundance", x = "Sample") +
  coord_cartesian(ylim = c(0,1), expand=0) + 
  scale_color_manual(values=SteppedSequential5Steps) + 
  scale_fill_manual(values=SteppedSequential5Steps) +
  facet_wrap(~Sample_type) + 
  theme(axis.text.x=element_blank())  -> p.fecal

p.fecal

ps_melt(ps.lt.phylum.ts) %>%
  ggplot(aes(fill=Phylum, y=Abundance, x=Sample_Order)) + 
    geom_bar(position="fill", stat="identity", width = 1) +
  theme_bw()  +   
  labs(y="Relative Abundance", x = "Sample") +
  coord_cartesian(ylim = c(0,1), expand=0) + 
  scale_color_manual(values=SteppedSequential5Steps) + 
  scale_fill_manual(values=SteppedSequential5Steps) +
  facet_wrap(~Sample_type) + 
  theme(axis.text.x=element_blank())  -> p.lt

p.lt

ragg::agg_jpeg(file.path(fig.dir, "fig_RelAb_fecal.jpeg"), width = 10, height = 6, units = "in", res = 600)
p.fecal
dev.off()
## png 
##   2
ragg::agg_jpeg(file.path(fig.dir, "fig_RelAb_lt.jpeg"), width = 5, height = 6, units = "in", res = 600)
p.lt
dev.off()
## png 
##   2

4 Alpha Diversity: LT Samples

4.1 Calculate

ps.LT.asv %>% estimate_richness(measures = c("Observed", "Shannon")) -> LT_richness
LT_richness$Faith <- pd(otu_table(ps.LT.asv), phy_tree(ps.LT.asv), include.root = FALSE)$PD 
LT_richness %>% rownames_to_column(var = "Label") -> LT_richness
gsub("\\.", "-", LT_richness$Label) -> LT_richness$Label 

meta.df %>%  
  filter(meta.df$Label %in% LT_richness$Label) %>% # merge with metadata
  select(c(Label:DNA_ng_ul_raw, "Group")) %>% 
  right_join(LT_richness, by = "Label") -> LT_richness

LT_richness  %>% pivot_longer(
  cols = contains(c("Observed", "Shannon", "Faith")),
  names_to = c("Measurement")
) -> LT_richness_stat # dataframe that works well for ggplot2() and anova_test()

LT_richness_stat$Mouse <- factor(LT_richness_stat$Mouse)
LT_richness_stat$Measurement <- factor(LT_richness_stat$Measurement)
LT_richness_stat$Cohort <- factor(LT_richness_stat$Cohort)
LT_richness %>% nrow()
## [1] 44
LT_richness_stat %>% nrow()
## [1] 132
44*3
## [1] 132

4.2 No surgery mice

We are most interested in intranasal treatment, diet, and genotype if possible. We also note that there were some mice whose lung tissues were collected on week 10 vs. week 13.

LT_richness %>% filter(Surgery == "None") -> lt.richness.nosurgery 
lt.richness.nosurgery %>%
  dplyr::count(Week, Genotype, Intranasal_Treatment) %>% arrange(n)
lt.richness.nosurgery %>%
  dplyr::count(Week, Intranasal_Treatment, Diet) %>% arrange(n)

Let’s graph:

LT_richness_stat %>% filter(Surgery == "None") -> lt.richness.nosurgery.long

lt.richness.nosurgery.long %>% 
  ggplot(aes(x = Intranasal_Treatment, y = value)) +
  geom_boxplot() + geom_point() + 
  facet_wrap(vars(Measurement), scales = "free") + theme_bw()

lt.richness.nosurgery.long %>% 
  ggplot(aes(x = Genotype, y = value)) +
  geom_boxplot() + geom_point() + 
  facet_wrap(vars(Measurement), scales = "free") + theme_bw()

lt.richness.nosurgery.long %>% 
  ggplot(aes(x = Diet, y = value)) +
  geom_boxplot() + geom_point() + 
  facet_wrap(vars(Measurement), scales = "free") + theme_bw()

lt.richness.nosurgery.long %>% 
  ggplot(aes(x = as.character(Week), y = value)) +
  geom_boxplot() + geom_point() + 
  facet_wrap(vars(Measurement), scales = "free") + theme_bw()

lt.richness.nosurgery.long %>% 
  ggplot(aes(x = Genotype, y = value)) +
  geom_boxplot() + geom_point() + 
  facet_grid(vars(Measurement), vars(Intranasal_Treatment), 
             scales = "free") + theme_bw()

lt.richness.nosurgery.long %>% 
  ggplot(aes(x = Diet, y = value)) +
  geom_boxplot() + geom_point() + 
  facet_grid(vars(Measurement), vars(Intranasal_Treatment),
             scales = "free") + theme_bw()

lt.richness.nosurgery.long %>% 
  ggplot(aes(x = as.character(Week), y = value)) +
  geom_boxplot() + geom_point() + 
  facet_grid(vars(Measurement), vars(Intranasal_Treatment),
             scales = "free") + theme_bw()

Graphically, it looks like intranasal treatment, week, and genotype may be of interest.

lt.richness.nosurgery %>% nrow()
## [1] 26
lt.richness.nosurgery %>%
  dplyr::count(Week, Intranasal_Treatment) %>% arrange(n)
lt.richness.nosurgery.long %>% 
  group_by(Week, Intranasal_Treatment, Measurement) %>%
  shapiro_test(value) %>%
  add_significance()

Mainly indistinguishable from normal distribution.

lt.richness.nosurgery.long %>% 
  group_by(Week, Intranasal_Treatment, Measurement) %>%
  identify_outliers(value) -> lt.nosurgery.outliers

lt.nosurgery.outliers

Outliers only in week 10 groups.

meta.sam %>% filter(Sample_type == "Lung Tissue") %>% 
  filter(Label %in% lt.nosurgery.outliers$Label) %>% 
  dplyr::select(Label, Sample_type:DNA_ng_ul_raw, Notes)

Interesting, maybe genotype or sex is playing a role here?

lt.richness.nosurgery.long %>% 
  ggqqplot(x = "value") + facet_grid(vars(Measurement), vars(Intranasal_Treatment), scales = "free")

lt.richness.nosurgery.long %>% 
  group_by(Measurement) %>% 
  anova_test(
    dv = value,
    between = c(Week, Intranasal_Treatment, Genotype),
  ) %>% get_anova_table() %>%
  adjust_pvalue(method = "holm") %>% 
  add_significance()

Nothing is significant, sadly. Let’s graph with intranasal treatment, as that had the greatest effect on host respiratory health metrics.

lt.richness.nosurgery.long %>%
  group_by(Measurement) %>%
  t_test(value ~ Intranasal_Treatment) %>% 
  adjust_pvalue(method = "holm") %>% 
  add_significance() -> lt.richness.nosurgery.res

lt.richness.nosurgery.res
# Change facet label names
alpha.labs <- c("Faith's PD", "Observed Richness", "Shannon")
names(alpha.labs) <- c("Faith", "Observed", "Shannon")

#lt.richness.nosurgery.res$y.position <- c(40, 225, 4.7)

lt.richness.nosurgery.long %>%
  ggplot(aes(x = Intranasal_Treatment, y = value, color = Intranasal_Treatment))+
  geom_boxplot() + 
  geom_point() + 
  facet_wrap(~Measurement, scales = "free_y",
             labeller = labeller(Measurement = alpha.labs)) + expand_limits(y = 0) + 
  theme_bw(base_size = 14) + scale_fill_brewer(palette="RdBu") +
  labs(y="Alpha Diversity Value", x="Intranasal Treatment", color = "Intranasal Treatment") + 
  theme(legend.position = "top") -> lt.alpha.group1.week
#  stat_pvalue_manual(lt.richness.nosurgery.res, label = "p.adj.signif")    -> lt.alpha.group1.week

lt.alpha.group1.week

4.3 HFD: effect of surgery?

For the after surgery timepoint for lung tissue, look at week 13 only and HFD.

LT_richness %>% filter(Diet == "High_Fat_Diet" & Week == 13) %>%
  dplyr::count(Surgery, Intranasal_Treatment)
LT_richness %>% filter(Diet == "High_Fat_Diet" & Week == 13) %>%
  dplyr::count(Surgery)
LT_richness %>% filter(Diet == "High_Fat_Diet" & Week == 13) %>%
  dplyr::count(Intranasal_Treatment)
LT_richness %>% filter(Diet == "High_Fat_Diet" & Week == 13)%>% nrow()
## [1] 23
LT_richness_stat %>% filter(Diet == "High_Fat_Diet" & Week == 13) -> lt.richness.hfd
lt.richness.hfd  %>%
  group_by(Measurement, Surgery) %>%
  shapiro_test(value) %>% 
  add_significance()
lt.richness.hfd %>% ggqqplot(x = "value") + facet_grid(vars(Measurement), vars(Surgery), scales = "free")

It looks like this is due to an outlier, but rarefaction curves looked like they plateaued. Proceed with caution.

lt.richness.hfd %>% 
  group_by(Measurement, Surgery) %>%
  shapiro_test(value) %>% 
  add_significance()
lt.richness.hfd %>% 
  group_by(Measurement, Surgery) %>%
  identify_outliers(value) 
lt.richness.hfd %>%
  ggplot(aes(x = Surgery, y = value)) + 
  geom_boxplot() + geom_point() + 
  facet_wrap(vars(Measurement), scales = "free") + theme_bw()

lt.richness.hfd %>% 
  group_by(Measurement) %>% 
  anova_test(
    dv = value,
    between = c(Surgery, Intranasal_Treatment),
  ) %>% get_anova_table() %>% 
  adjust_pvalue(method = "holm") %>% 
  add_significance()

Nothing is significant. Let’s make a graph showing surgery:

lt.richness.hfd %>%
  ggplot(aes(x = Surgery, y = value, color = Surgery)) + 
  geom_boxplot() + geom_point() + 
  facet_wrap(~Measurement, scales = "free_y",
             labeller = labeller(Measurement = alpha.labs)) + expand_limits(y = 0) + 
  theme_bw(base_size = 14) + scale_color_brewer(palette="Dark2") +
  labs(y="Alpha Diversity Value", x="Surgery", color = "Surgery") + 
  theme(legend.position = "top") -> fig.lt.hfd

fig.lt.hfd

4.4 Summary figure

plot_grid(
  lt.alpha.group1.week, fig.lt.hfd,
  labels = "AUTO", ncol = 2, label_size = 20, scale = 0.95,
  rel_widths = c(1, 1.3)
) -> fig.lt.alpha

fig.lt.alpha

ragg::agg_jpeg(file.path(fig.dir, "fig_lt_alpha.jpeg"), width = 14, height = 5, units = "in", res = 600, scaling = 1)
fig.lt.alpha
dev.off()
## png 
##   2

5 Lung tissue beta diversity

5.1 No surgery mice

For PERMANOVA, let’s have n>=5 for each group.

ps.LT.asv.ts %>% subset_samples(Surgery == "None")  -> ps.lt.ts.nosurgery

# Calculate ditances
ps.lt.ts.nosurgery %>% phyloseq::distance(method = "unifrac") -> LT_nosurgery_uni
ps.lt.ts.nosurgery %>% phyloseq::distance(method = "wunifrac") -> LT_nosurgery_weighted_uni

# Ordinate distances
LT_nosurgery.uni_ord <- ordinate(ps.lt.ts.nosurgery, method = "PCoA", distance = "unifrac")
LT_nosurgery.weighted_uni_ord <- ordinate(ps.lt.ts.nosurgery, method = "PCoA", distance = "wunifrac")

# Extract dataframe
ps.lt.ts.nosurgery@sam_data %>% data.frame() -> LT_nosurgery.df

LT_nosurgery.df %>% dplyr::count(Genotype, Intranasal_Treatment) %>% arrange(n)
# Set seed for reproducibility
# UniFrac
set.seed(123)

adonis2(formula = LT_nosurgery_uni ~ Intranasal_Treatment * Genotype, data = LT_nosurgery.df, permutations = 9999, parallel = getOption("mc.cores"))
# Weighted UniFrac
set.seed(123)
adonis2(formula = LT_nosurgery_weighted_uni ~ Intranasal_Treatment  * Genotype, data = LT_nosurgery.df, permutations = 9999, parallel = getOption("mc.cores")) 

Not significant.

5.1.1 Graphs

plot_ordination(ps.lt.ts.nosurgery, LT_nosurgery.uni_ord, type = "samples", color = "Intranasal_Treatment") + theme_bw()  + 
  labs(title = "UniFrac") +
  scale_color_brewer(palette="Dark2")

plot_ordination(ps.lt.ts.nosurgery, LT_nosurgery.weighted_uni_ord, type = "samples", color = "Intranasal_Treatment") + theme_bw()  + 
  labs(title = "Weighted UniFrac") +
  scale_color_brewer(palette="Dark2")

plot_ordination(ps.lt.ts.nosurgery, LT_nosurgery.uni_ord, type = "samples", color = "Intranasal_Treatment") + theme_bw(base_size = 14)  + 
  scale_color_brewer(palette="Dark2") +
  theme(legend.position = "top") -> p.lt.nosurgery.uni
p.lt.nosurgery.uni

plot_ordination(ps.lt.ts.nosurgery, LT_nosurgery.weighted_uni_ord, type = "samples", color = "Intranasal_Treatment") + 
  theme_bw(base_size = 14)  + 
  scale_color_brewer(palette="Dark2") +
  theme(legend.position = "top") -> p.lt.nosurgery.wuni
p.lt.nosurgery.wuni

plot_grid(p.lt.nosurgery.uni, p.lt.nosurgery.wuni,
          ncol = 2, labels = "AUTO", label_size = 20) -> p.lt.nosurgery

p.lt.nosurgery

ragg::agg_jpeg(file.path(fig.dir, "fig_lt_nosurgery_beta.jpeg"), width = 10, height = 5, units = "in", res = 600, scaling = 1.1)
p.lt.nosurgery
dev.off()
## png 
##   2

5.2 Did Surgery result in lung microbiome changes?

Only look at at HFD at week 13.

ps.LT.asv.ts %>% subset_samples(Diet == "High_Fat_Diet") %>% subset_samples(Week == "13")  -> ps.lt.ts.hfd

# Calculate ditances
ps.lt.ts.hfd %>% phyloseq::distance(method = "unifrac") -> LT_hfd_uni
ps.lt.ts.hfd %>% phyloseq::distance(method = "wunifrac") -> LT_hfd_weighted_uni

# Ordinate distances
LT_hfd.uni_ord <- ordinate(ps.lt.ts.hfd, method = "PCoA", distance = "unifrac")
LT_hfd.weighted_uni_ord <- ordinate(ps.lt.ts.hfd, method = "PCoA", distance = "wunifrac")

# Extract dataframe
ps.lt.ts.hfd@sam_data %>% data.frame() -> LT_hfd.df

LT_hfd.df %>% dplyr::count(Surgery) %>% arrange(n)
# Set seed for reproducibility
# UniFrac
set.seed(123)

adonis2(formula = LT_hfd_uni ~ Surgery, data = LT_hfd.df, permutations = 9999, parallel = getOption("mc.cores"))
# Weighted UniFrac
set.seed(123)
adonis2(formula = LT_hfd_weighted_uni ~ Surgery, data = LT_hfd.df, permutations = 9999, parallel = getOption("mc.cores")) 

NS.

plot_ordination(ps.lt.ts.hfd, LT_hfd.uni_ord, type = "samples", color = "Surgery") + theme_bw(base_size = 14)  + 
  scale_color_brewer(palette="Dark2") +
  theme(legend.position = "top") -> p.lt.hfd.uni
p.lt.hfd.uni

plot_ordination(ps.lt.ts.hfd, LT_hfd.weighted_uni_ord, type = "samples", color = "Surgery") + 
  theme_bw(base_size = 14)  + 
  scale_color_brewer(palette="Dark2") +
  theme(legend.position = "top") -> p.lt.hfd.wuni
p.lt.hfd.wuni

plot_grid(p.lt.hfd.uni, p.lt.hfd.wuni,
          ncol = 2, labels = "AUTO", label_size = 20) -> p.lt.hfd

p.lt.hfd

ragg::agg_jpeg(file.path(fig.dir, "fig_lt_hfd_beta.jpeg"), width = 10, height = 5, units = "in", res = 600, scaling = 1.1)
p.lt.hfd
dev.off()
## png 
##   2

6 Alpha Diversity: Fecal Samples

6.1 Calculate

ps.F.asv
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3304 taxa and 177 samples ]
## sample_data() Sample Data:       [ 177 samples by 40 sample variables ]
## tax_table()   Taxonomy Table:    [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq()      DNAStringSet:      [ 3304 reference sequences ]
ps.F.asv %>% estimate_richness(measures = c("Observed", "Shannon")) -> F_richness
## Warning in estimate_richness(., measures = c("Observed", "Shannon")): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
## 
## We recommended that you find the un-trimmed data and retry.
F_richness$Faith <- pd(otu_table(ps.F.asv), phy_tree(ps.F.asv), include.root = FALSE)$PD 
F_richness %>% rownames_to_column(var = "Label") -> F_richness
gsub("X", "", F_richness$Label) -> F_richness$Label # remove X in front of sample names

meta.df %>%  
  filter(meta.df$Label %in% F_richness$Label) %>% # merge with metadata
  select(c(Label, Genotype:Surgery, Kit:Feces_weight, Week, Cohort, Mouse, Notes)) %>% 
  right_join(F_richness, by = "Label") -> F_richness

F_richness %>% group_by(Mouse, Week) %>% filter(n()>1) -> F_richness.dup

F_richness.dup
# calculate averages for duplicates
F_richness.dup %>%
  dplyr::select(Label, Mouse, Week, Observed:Faith) %>% group_by(Mouse) %>%
  dplyr::summarise(Observed = mean(Observed),
                   Shannon = mean(Shannon),
                   Faith = mean(Faith)) -> F_richness.dup.calc

# add week information for duplicates
F_richness.dup %>% dplyr::select(Mouse, Week) %>% unique() %>%
  ungroup() %>%
  right_join(F_richness.dup.calc, by = "Mouse") -> F_richness.dup.calc 

# add metadata for duplicates
meta.df %>% filter(Sample_type == "Fecal") %>%
  filter(Mouse %in% F_richness.dup.calc$Mouse) %>% 
  select(c(Genotype:Surgery, Week, Cohort, Mouse, Notes)) %>% 
  right_join(F_richness.dup.calc, by = c("Mouse", "Week")) %>%
  unique() -> F_richness.dup.calc

F_richness.dup.calc
# Clean up: remove duplicates from original dataframe, then add calculated values. 
F_richness %>% filter(!Label %in% F_richness.dup$Label) %>%
  bind_rows(F_richness.dup.calc) -> F_richness
F_richness %>% filter(is.na(Diet) == FALSE) -> F_richness

F_richness  %>% pivot_longer(
  cols = contains(c("Observed", "Shannon", "Faith")),
  names_to = c("Measurement")
) -> F_richness_stat # dataframe that works well for ggplot2() and anova_test()

F_richness_stat$Mouse <- factor(F_richness_stat$Mouse)
F_richness_stat$Measurement <- factor(F_richness_stat$Measurement)
F_richness_stat$Cohort <- factor(F_richness_stat$Cohort)

6.2 No Surgery group: effect of diet?

I want to see if high fat diet resulted in a significant decline in alpha diversity, as has been reported before. To do this, I will compare wk 0 and wk 10 alpha diversity for all mice based on diet (to increase n). I also want to see if genotype and sex had an effect.

F_richness %>% filter(Surgery == "None") -> F_richness.nosurgery
F_richness.nosurgery %>%
  dplyr::count(Week, Diet, Genotype, Sex) %>% arrange(n)
F_richness.nosurgery %>%
  dplyr::count(Week, Diet, Genotype) %>% arrange(n)
F_richness.nosurgery %>% distinct(Mouse) %>% nrow()
## [1] 58
F_richness.nosurgery %>% nrow()
## [1] 117
F_richness_stat %>% filter(Surgery == "None") -> F_richness.nosurgery.long

F_richness.nosurgery.long %>%
  group_by(Measurement, Week, Diet, Genotype) %>% 
  shapiro_test(value) %>%
  add_significance() -> F_richness.nosurgery.shapirowilk

F_richness.nosurgery.shapirowilk
F_richness.nosurgery.shapirowilk %>% filter(p<0.05)

All indistinguishable from normal distribution.

F_richness.nosurgery.long %>%
  group_by(Measurement, Week, Diet, Genotype) %>% 
  identify_outliers(value)

The feces weight, gDNA concentrations, and notes (or the lack thereof) for even extreme outliers look good.

F_richness.nosurgery.long %>%
  group_by(Measurement) %>% 
  anova_test(dv = value, 
             wid = Mouse, 
             within = Week,
             between = c(Diet, Genotype)) %>%
  get_anova_table() %>% 
  as.data.frame() %>% 
  adjust_pvalue(method = "holm") %>% 
  add_significance() -> F.nosurgery.res

F.nosurgery.res
F.nosurgery.res %>% filter(p.adj < 0.05)

Let’s graph what they look like:

# Change facet label names
alpha.labs <- c("Faith's PD", "Observed Richness", "Shannon")
names(alpha.labs) <- c("Faith", "Observed", "Shannon")

ggplot(F_richness.nosurgery.long, aes(x = as.factor(Week), y = value, color = Diet))+
  geom_boxplot(position=position_dodge(width = 0.8)) + geom_point(position=position_dodge(width = 0.8)) + 
  facet_wrap(~Measurement, scales = "free", labeller = labeller(Measurement = alpha.labs)) + 
  labs(x="Week", y="Alpha Diversity Metric Value", fill="Diet")  + expand_limits(y = 0) + theme_bw() +
  theme(legend.position = "top") 

# plot mean +/- standard error
fun_se <- function(x){
  c(std_error = sd(x)/sqrt(length(x)), mean = mean(x))
}

# calculate mean and se for each week, measurement, and diet
F_richness.nosurgery.long %>%
  doBy::summary_by(value ~ Measurement + Week + Diet, FUN = fun_se) -> F_richness.nosurgery.graph

ggplot(F_richness.nosurgery.graph, aes(x = Week, y = value.mean, color = Diet))+
  geom_line(aes(group=Diet)) + 
  geom_errorbar(aes(ymin=value.mean-value.std_error, ymax=value.mean+value.std_error, width=1)) +
  theme_bw(base_size = 14) + scale_fill_brewer(palette="RdBu") + 
  facet_wrap(~Measurement, scales = "free_y", 
             labeller = labeller(Measurement = alpha.labs)) + 
  labs(x="Week", y="Alpha Diversity Metric Value", fill="Diet")  + expand_limits(y = 0) +
  theme(legend.position = "top") -> fig.alpha.nosurgery

fig.alpha.nosurgery

ragg::agg_jpeg(file.path(fig.dir, "fig.alpha.nosurgery.jpeg"), width = 6, height = 5, units = "in", res = 600)
fig.alpha.nosurgery
dev.off()
## png 
##   2

6.3 Effect of surgery?

F_richness %>%
  filter(Diet == "High_Fat_Diet" & Week == 13) -> f.richness.hfd  

f.richness.hfd %>% dplyr::count(Surgery) %>% arrange(n)
F_richness %>%
  filter(Diet == "High_Fat_Diet" & Week == 13)  %>% distinct(Mouse) %>% nrow()
## [1] 45
f.richness.hfd %>% dplyr::count(Surgery, Genotype) %>% arrange(n)
f.richness.hfd %>% dplyr::count(Surgery, Genotype, Sex) %>% arrange(n)
F_richness_stat %>%
  filter(Diet == "High_Fat_Diet" & Week == 13) -> f.richness.hfd
f.richness.hfd %>%
  group_by(Measurement, Surgery, Genotype) %>% 
  shapiro_test(value) %>%
  add_significance()
f.richness.hfd %>%
  group_by(Measurement, Surgery, Genotype) %>% 
  identify_outliers(value)
f.richness.hfd %>% filter(Genotype == "KO") %>% 
  ggqqplot(x = "value") + facet_grid(vars(Measurement), vars(Surgery), scales = "free")

The QQ plots look quite good, except for the few outliers.

f.richness.hfd %>%
  group_by(Measurement) %>% 
  anova_test(dv = value, 
             wid = Mouse, 
             between = c(Surgery, Genotype, Sex)) %>%
  get_anova_table() %>% 
  as.data.frame() %>% 
  adjust_pvalue(method = "holm") %>% 
  add_significance() -> F.hfd.res

F.hfd.res

All ns. Let’s graph:

f.richness.hfd %>%
  ggplot(aes(x=Surgery, y=value, color=Surgery)) + 
  geom_boxplot() + 
  geom_point() + 
  scale_fill_brewer(palette="RdBu") + 
  labs( y="Alpha Diversity Metric value", color="Surgery", x = "Surgery") + 
  facet_wrap(vars(Measurement), scales = "free", labeller = labeller(Measurement = alpha.labs)) +
  theme_bw(base_size = 14)  + expand_limits(y = 0) + theme(legend.position = "top") -> fig.alpha.hfd
fig.alpha.hfd

ragg::agg_jpeg(file.path(fig.dir, "fig_f_alpha_hfd.jpeg"), width = 7, height = 5, units = "in", res = 600)
fig.alpha.hfd
dev.off()
## png 
##   2

7 Beta Diversity with fecal samples: Diet

ps.F.asv.ts %>% subset_samples(Surgery == "None") -> F_Diet

F_Diet %>%
  phyloseq::distance(method = "unifrac") -> F_uni
F_Diet %>%
  phyloseq::distance(method = "wunifrac") -> F_weighted_uni

F_Diet@sam_data %>% data.frame() -> F.df

F_uni_ord <- ordinate(F_Diet, method = "PCoA", distance = "unifrac")
F_weighted_uni_ord <- ordinate(F_Diet, method = "PCoA", distance = "wunifrac")
F.df %>% dplyr::count(Diet, Week) %>% arrange(n)
F.df %>% dplyr::count(Diet, Genotype, Week) %>% arrange(n)
F.df %>% nrow()
## [1] 120
F.df %>% distinct(Mouse) %>%  nrow()
## [1] 58
# Set seed for reproducibility
set.seed(123)

adonis2(formula = F_uni ~ Diet * Genotype * Week, data = F.df, permutations = 9999, parallel = getOption("mc.cores")) -> F_diet_uni.res
F_diet_uni.res
set.seed(123)
adonis2(formula = F_weighted_uni ~ Diet * Genotype * Week, data = F.df, permutations = 9999, parallel = getOption("mc.cores")) -> F_diet_weighted_uni.res

F_diet_weighted_uni.res

As expected, diet is highly significant. Genotype is significant but not any of its interaction terms.

F.df %>% dplyr::count(Genotype, Diet) %>% arrange(n)
plot_ordination(F_Diet, F_uni_ord, type = "samples", color = "Diet") + theme_bw(base_size = 14)   + 
  scale_fill_brewer(palette="RdBu") + theme(legend.position = "top") + labs(title = "UniFrac distance") + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  facet_wrap(~Week) + stat_ellipse() -> fig.diet.uni

plot_ordination(F_Diet, F_weighted_uni_ord, type = "samples", color = "Diet") + 
  theme_bw(base_size = 14)   + 
  scale_fill_brewer(palette="RdBu") + theme(legend.position = "top", panel.spacing.x = unit(10, "mm")) + 
  labs(title = "Weighted UniFrac distance") + theme(plot.title = element_text(hjust = 0.5)) + 
  facet_wrap(~Week) + stat_ellipse() -> fig.diet.w.uni

fig.diet.uni 

fig.diet.w.uni

plot_ordination(F_Diet, F_uni_ord, type = "samples", color = "Genotype") + theme_bw(base_size = 14)   + labs(title = "UniFrac distance") + 
  scale_fill_brewer(palette="RdBu") + theme(legend.position = "top", plot.title = element_text(hjust = 0.5))  + stat_ellipse() -> fig.geno.uni

plot_ordination(F_Diet, F_weighted_uni_ord, type = "samples", color = "Genotype") + labs(title = "Weighted UniFrac distance") + 
  theme_bw(base_size = 14)   + 
  scale_fill_brewer(palette="Dark2") + theme(legend.position = "top", panel.spacing.x = unit(8, "mm"), plot.title = element_text(hjust = 0.5)) + stat_ellipse() -> fig.geno.w.uni

# UniFrac: diet*genotype was significant
plot_ordination(F_Diet, F_uni_ord, type = "samples", color = "Genotype") + theme_bw(base_size = 14)   + 
  scale_color_brewer(palette="Dark2") + facet_wrap(~Diet) +  labs(title = "UniFrac distance") + 
  theme(legend.position = "top", plot.title = element_text(hjust = 0.5))  + stat_ellipse() -> fig.diet.geno.uni

fig.geno.uni

fig.diet.geno.uni

fig.geno.w.uni

plot_ordination(F_Diet, F_uni_ord, type = "samples", color = "Genotype") + theme_bw(base_size = 14)   + 
  scale_fill_brewer(palette="RdBu") +
  theme(legend.position = "top", panel.spacing.x = unit(5, "mm")) +
  facet_grid(vars(Diet),vars(Week)) + stat_ellipse()  -> fig.diet.geno.week.uni
fig.diet.geno.week.uni

plot_grid(
  fig.alpha.nosurgery,fig.diet.uni,
  fig.diet.w.uni, fig.diet.geno.uni, labels = "AUTO", ncol = 2, label_size = 20
) -> fig.diet.summary

fig.diet.summary

ragg::agg_jpeg(file.path(fig.dir, "fig_diet_summary.jpeg"), width = 18, height = 12, units = "in", res = 600)
fig.diet.summary
dev.off()
## png 
##   2
plot_grid(
  fig.geno.uni, fig.geno.w.uni, labels = "AUTO", ncol = 2, label_size = 20
) -> fig.nosurgery.geno.summary

fig.nosurgery.geno.summary

ragg::agg_jpeg(file.path(fig.dir, "fig_geno_summary.jpeg"), width = 12, height = 6, units = "in", res = 600)
fig.nosurgery.geno.summary
dev.off()
## png 
##   2

7.1 Fecal beta diversity: HFD

HFD: Effect of surgery; focus on wk 13 only.

ps.F.asv.ts %>% subset_samples(Diet == "High_Fat_Diet" & Week == 13) -> F_HFD
F_HFD
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3304 taxa and 45 samples ]
## sample_data() Sample Data:       [ 45 samples by 40 sample variables ]
## tax_table()   Taxonomy Table:    [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq()      DNAStringSet:      [ 3304 reference sequences ]
F_HFD %>%
  phyloseq::distance(method = "unifrac") -> F_uni.HFD
F_HFD %>%
  phyloseq::distance(method = "wunifrac") -> F_weighted_uni.HFD

F_HFD@sam_data %>% data.frame() -> F_HFD.df

F_uni_ord.HFD <- ordinate(F_HFD, method = "PCoA", distance = "unifrac")
F_weighted_uni_ord.HFD <- ordinate(F_HFD, method = "PCoA", distance = "wunifrac")
F_HFD.df %>% dplyr::count(Surgery) %>% arrange(n)
F_HFD.df %>% dplyr::count(Genotype) %>% arrange(n)
F_HFD.df %>% dplyr::count(Surgery, Genotype) %>% arrange(n)
F_HFD.df %>% nrow()
## [1] 45
# Set seed for reproducibility
set.seed(123)

adonis2(formula = F_uni.HFD ~ Surgery * Genotype, data = F_HFD.df, permutations = 9999, parallel = getOption("mc.cores")) 
set.seed(123)
adonis2(formula = F_weighted_uni.HFD ~ Surgery * Genotype, data = F_HFD.df, permutations = 9999, parallel = getOption("mc.cores")) 

Genotype is significant, surgery is not.

plot_ordination(F_Diet, F_weighted_uni_ord, type = “samples”, color = “Genotype”) + labs(title = “Weighted UniFrac distance”) + theme_bw(base_size = 14) + scale_fill_brewer(palette=“RdBu”) + theme(legend.position = “top”, panel.spacing.x = unit(8, “mm”), plot.title = element_text(hjust = 0.5)) + stat_ellipse() -> fig.geno.w.uni

# stat_ellipse for the ones that are * 
plot_ordination(F_HFD, F_uni_ord.HFD, type = "samples", color = "Genotype") + theme_bw(base_size = 14)   + 
  labs(title = "UniFrac distance") + 
  theme(legend.position = "top", plot.title = element_text(hjust = 0.5)) +  stat_ellipse() -> fig.hfd.uni.geno
fig.hfd.uni.geno

plot_ordination(F_HFD, F_weighted_uni_ord.HFD, type = "samples", color = "Genotype") + theme_bw(base_size = 14)   + 
  labs(title = "Weighted UniFrac distance") + 
  theme(legend.position = "top", plot.title = element_text(hjust = 0.5)) + 
  stat_ellipse() -> fig.hfd.w.uni.geno
fig.hfd.w.uni.geno

# Graphing out of curiosity  
plot_ordination(F_HFD, F_uni_ord.HFD, type = "samples", color = "Surgery") + theme_bw(base_size = 14)   + 
  scale_color_brewer(palette="Dark2") + theme(legend.position = "top") + stat_ellipse()

plot_ordination(F_HFD, F_weighted_uni_ord.HFD, type = "samples", color = "Surgery") + theme_bw(base_size = 14)   + 
  scale_color_brewer(palette="Dark2") + theme(legend.position = "top")  + stat_ellipse()

plot_grid(
  fig.hfd.uni.geno,
  fig.hfd.w.uni.geno, 
  labels = "AUTO", ncol = 2, label_size = 20
) -> fig.hfd.summary.beta

fig.hfd.summary.beta

ragg::agg_jpeg(file.path(fig.dir, "fig_hfd_beta.jpeg"), width = 12, height = 5, units = "in", res = 600)
fig.hfd.summary.beta
dev.off()
## png 
##   2

8 Reproducibility

sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] ggtext_0.1.2                cowplot_1.1.1              
##  [3] ragg_1.2.6                  doBy_4.6.19                
##  [5] DESeq2_1.40.2               SummarizedExperiment_1.30.2
##  [7] Biobase_2.60.0              MatrixGenerics_1.12.3      
##  [9] matrixStats_1.1.0           GenomicRanges_1.52.1       
## [11] GenomeInfoDb_1.36.4         IRanges_2.34.1             
## [13] S4Vectors_0.38.2            BiocGenerics_0.46.0        
## [15] phytools_1.9-16             maps_3.4.1                 
## [17] picante_1.8.2               nlme_3.1-163               
## [19] ape_5.7-1                   plotly_4.10.3              
## [21] DT_0.29                     ggrepel_0.9.4              
## [23] microbiome_1.22.0           colorBlindness_0.1.9       
## [25] RColorBrewer_1.1-3          pals_1.8                   
## [27] microViz_0.11.0             ggpubr_0.6.0               
## [29] rstatix_0.7.2               vegan_2.6-4                
## [31] lattice_0.22-5              permute_0.9-7              
## [33] lubridate_1.9.3             forcats_1.0.0              
## [35] stringr_1.5.0               dplyr_1.1.3                
## [37] purrr_1.0.2                 readr_2.1.4                
## [39] tidyr_1.3.0                 tibble_3.2.1               
## [41] ggplot2_3.4.4               tidyverse_2.0.0            
## [43] fs_1.6.3                    phyloseq_1.44.0            
## [45] magrittr_2.0.3             
## 
## loaded via a namespace (and not attached):
##   [1] rstudioapi_0.15.0       jsonlite_1.8.7          farver_2.1.1           
##   [4] rmarkdown_2.25          zlibbioc_1.46.0         vctrs_0.6.4            
##   [7] multtest_2.56.0         RCurl_1.98-1.13         S4Arrays_1.2.1         
##  [10] htmltools_0.5.6.1       plotrix_3.8-2           broom_1.0.5            
##  [13] Rhdf5lib_1.22.1         rhdf5_2.44.0            gridGraphics_0.5-1     
##  [16] sass_0.4.7              bslib_0.5.1             htmlwidgets_1.6.2      
##  [19] plyr_1.8.9              cachem_1.0.8            igraph_1.5.1           
##  [22] lifecycle_1.0.3         iterators_1.0.14        pkgconfig_2.0.3        
##  [25] Matrix_1.6-1.1          R6_2.5.1                fastmap_1.1.1          
##  [28] GenomeInfoDbData_1.2.10 digest_0.6.33           numDeriv_2016.8-1.1    
##  [31] colorspace_2.1-0        crosstalk_1.2.0         textshaping_0.3.7      
##  [34] labeling_0.4.3          clusterGeneration_1.3.8 fansi_1.0.5            
##  [37] timechange_0.2.0        httr_1.4.7              abind_1.4-5            
##  [40] mgcv_1.9-0              compiler_4.3.1          microbenchmark_1.4.10  
##  [43] bit64_4.0.5             withr_2.5.1             doParallel_1.0.17      
##  [46] backports_1.4.1         BiocParallel_1.34.2     optimParallel_1.0-2    
##  [49] carData_3.0-5           ggsignif_0.6.4          MASS_7.3-60            
##  [52] DelayedArray_0.26.7     scatterplot3d_0.3-44    biomformat_1.28.0      
##  [55] tools_4.3.1             glue_1.6.2              quadprog_1.5-8         
##  [58] rhdf5filters_1.12.1     gridtext_0.1.5          grid_4.3.1             
##  [61] Rtsne_0.16              cluster_2.1.4           reshape2_1.4.4         
##  [64] ade4_1.7-22             generics_0.1.3          gtable_0.3.4           
##  [67] tzdb_0.4.0              data.table_1.14.8       hms_1.1.3              
##  [70] xml2_1.3.5              Deriv_4.1.3             car_3.1-2              
##  [73] utf8_1.2.4              XVector_0.40.0          foreach_1.5.2          
##  [76] pillar_1.9.0            vroom_1.6.4             splines_4.3.1          
##  [79] bit_4.0.5               survival_3.5-7          tidyselect_1.2.0       
##  [82] locfit_1.5-9.8          Biostrings_2.68.1       knitr_1.44             
##  [85] xfun_0.40               expm_0.999-7            stringi_1.7.12         
##  [88] lazyeval_0.2.2          yaml_2.3.7              evaluate_0.22          
##  [91] codetools_0.2-19        cli_3.6.1               systemfonts_1.0.5      
##  [94] munsell_0.5.0           jquerylib_0.1.4         dichromat_2.0-0.1      
##  [97] Rcpp_1.0.11             mapproj_1.2.11          coda_0.19-4            
## [100] ellipsis_0.3.2          bitops_1.0-7            phangorn_2.11.1        
## [103] viridisLite_0.4.2       scales_1.2.1            crayon_1.5.2           
## [106] combinat_0.0-8          rlang_1.1.1             fastmatch_1.1-4        
## [109] mnormt_2.1.1